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0^ ' A given set of data-points in some feature space may be associated with a Schrodinger equation 

, whose potential is determined by the data. This is known to lead to good clustering solutions. Here 

■ we extend this approach into a full-fledged dynamical scheme using a time-dependent Schrodinger 

(N ■ equation. Moreover, we approximate this Hamiltonian formalism by a truncated calculation within 

a set of Gaussian wave functions (coherent states) centered around the original points. This al- 
1-^ ■ lows for analytic evaluation of the time evolution of all such states, opening up the possibility of 

exploration of relationships among data-points through observation of varying dynamical-distances 
among points and convergence of points into clusters. This formalism may be further supplemented 
00 ■ by preprocessing, such as dimensional reduction through singular value decomposition or feature 

T-H ' filtering. 

PACS numbers: 89.75.Fb,89.90.+n,95.75.Pq,89.75.Kd 

? ' 

^ ! INTRODUCTION 

^ : 

Clustering of data is, in general, an ill-defined problem. Nonetheless it is a very important one in many scientific 
and technological fields of study. Given a set of data-points one looks for possible structures by sorting out which 
points are close to each other and, therefore, in some sense belong together. This is a preliminary stage taken before 



< 



investigating what properties are common to these subsets of the data. 

It has recently become quite popular to investigate such questions in a dynamic framework, thus allowing for more 
Q-i fluid associations of data-points, rather than rigid clusters. Diffusion geometry 0, 0) Q is such a method, based on a 
discrete analog of the heat equation 

^ ■ 5$ 

^ ; . 

, where H is some operator with positive eigenvalues, guaranteeing that the temporal evolution of t) is that of 
diffusion. Thus, starting out with <I>(x, 0), e.g. a Gaussian concentrated around some data point one would expect 
, t) to spread over all space that is occupied by the data points. 

' In this paper we advocate the use of a Schrodinger Hamiltonian H that is intimately connected to the data- 
structure, as defined by the quantum clustering method 5] and summarized below. We extend it into a time-dependent 
Schrodinger equation: 

.t: ^'-^=H^i.^t) (2) 

5_j ■ The ensuing Dynamic Quantum Clustering (DQC) formalism allows us, by varying a few parameters, to study in 
^ detail the temporal evolution of wave-functions representing the original data-points. In turn, this dynamical behavior 
allows us to explore the structure of the quantum potential function defined by the quantum clustering method. 

DQC begins by associating each data-point with a state in Hilbert space. The temporal development of the centroids 
of these states may be viewed in the original data-space as moving images of the original points. Their distances 
to one-another change with time, thus representing associations they form with each other. Convergence of many 
points onto a common center at some instant of time is an obvious manifestation of clustering. Many transitional 
relationships may occur, revealing substructures in clusters or even more complex associations. For this reason we 
propose this approach as a general method for visually and interactively searching for and exploring structures in sets 
of data. 



* This work was supported by the U. S. DOE, Contract No. DE-AC02-76SF00515. 
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Quantum Clustering 

The quantum clustering approach begins, as does the well-known Parzen- window estimator [3], by associating to 

each of n data points Xi in a Euclidean space of d dimensions a Gaussian wave-function -0^ (x) = e and then 

constructing the sum of all these Gaussians, 

,5-^ jX-Xj)'^ 

i,{x)=Y.^ (3) 

i 

Conventional scale-space clustering [llj views this function as a probability distribution (up to an overall factor) that 
could have generated the observed points, and regards therefore its maxima as determining locations of cluster centers. 
Often these maxima are not very prominent and, in order to uncover more of them, one has to reduce cr down to low 
values where the number and location of the maxima depend sensitively upon the choice of cr. 

Quantum clustering took a different approach, requiring to be the ground-state of the Hamiltonian 

2 

i/^ = (-yV2 + y(x))v = £;o^. (4) 

By positing this requirement |^], the potential function V^(x) has become inextricably bound to the system of data- 
points, since V{x) is determined, up to a constant, by a simple algebraic inversion of EqHl [l^ Moreover, one may 
expect V to have minima in regions where ip has maxima and furthermore, that these minima will be more pronounced 
than the corresponding maxima found in the Parzen estimator. In fact, it frequently turns out that a concentration 
of data-points will lead to a local minimum in V, even if does not display a local maximum. Thus, by replacing the 
problem of finding maxima of the Parzen estimator by the problem of locating the minima of the associated potential, 
V{x)^ we simplify the process of identifying clusters. The effectiveness of this approach has been demonstrated in 
the work by Horn and Gottlieb P, [13] ■ It should be noted that the enhancement of features obtained by applying 
Eq|3] comes from the interplay of two effects: attraction of the wave-function to the minima of V and spreading of 
the wave- function due to the second derivative (kinetic term). This may be viewed as an alternative model to the 
conventional probabilistic approach, incorporating attraction to cluster-centers and creation of noise, both inferred 
from - or realized by - the given experimental data. 

DQC drops the probabilistic interpretation of ij} and replaces it by that of a probability-amplitude, as customary in 
Quantum Mechanics. DQC is set up to associate data-points with cluster centers in a natural fashion. Whereas in QC 
this association was done by finding their loci on the slopes of V , here we follow the quantum-mechanical temporal 
evolvement of states associated with these points. Specifically, we will view each data-point as the expectation value 

(X-Xjf 

of the position operator in a Gaussian wave- function tpi{x) = e ^ ; the temporal development of this state traces 
the association of the data-point it represents with the minima of V{x) and thus, with the other data-points. 



Dynamic Quantum Clustering (DQC) 

As we already noted, the conversion of the static QC method to a full dynamical one, begins by focusing attention 

on the Gaussian wave- function, tl^iix) = C e 2^ , associated with the i data point, where C is the appropriate 
normalization factor. Thus, by construction, the expectation value of the operator x in this state is simply the 
coordinates of the original data point; i.e., 

X, = {ip\x\il;) = J dx^*{x) x^pi{x). (5) 

The dynamical part of the DQC algorithm is that, having constructed the potential function V{x), we study the time 
evolution of each state ipii^) ^ determined by the time dependent Schrodinger equation; i.e., 

^ H^,ix,t) ^ ( + V{x)) Ux,t), (6) 



dt ^ ' ' ' \ 2m 

where m is an arbitrarily chosen mass for a particle moving in (i-dimensions If we set m = then, by construction, 
i}{x) of Eq. 3 is the lowest energy eigenstate of the Hamiltonian. If m is chosen to have a different value, then not 
only does each individual state ^i{x) evolve in time, but so does the sum of the states, 'ip{x). 
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The important feature of quantum dynamics, which makes the evolution so useful in the clustering problem, is that 
according to Ehrenfest's theorem, the time-dependent expectation value 

{iP{t)\^\m)= J dxi;*{x,t)xM^,t), (7) 

satisfies the equation, 

d^{x{t)) 1 



= m)\^v{x)m))- (9) 

If Tpii^) is a narrow Gaussian, this is equivalent to saying that the center of each wave- function rolls towards the 
nearest minimum of the potential according to the classical Newton's law of motion. This means we can explore the 
relation of this data point to the minima of V{x) by following the time-dependent trajectory ( Xi{t) ) = {ipi{t)\ x\il)i{t)). 
Clearly, given Ehrenfest's theorem, we expect to sec any points located in, or near, the same local minimum of V{x) to 
oscillate about that minimum, coming together and moving apart. In our numerical solutions we generate animations 
which display this dynamics for a finite time. This allows us to visually trace the clustering of points associated with 
each one of the potential minima. 

In their quantum clustering paper Horn and Gottlieb successfully used classical gradient descent to cluster data by 
moving points (on classical trajectories) to the nearest local minimum of V{x). The idea being that points which end 
up at the same minimum are in the same cluster. At first glance it would seem that DQC replaces the conceptually 
simple problem of implementing gradient descent with the more difficult one of solving complicated partial differential 
equations. We will show the difficulty is only apparent. In fact, the solution of the Schrodinger equation can be 
simplified considerably, and will also allow further insights than the gradient descent method. The DQC algorithm 
translates the problem of solving the Schrodinger equation into a matrix form which captures most of the details of 
the analytic problem, but which involves N x A''-matrices whose dimension, N, is less than or equal to the number of 
data points. This reduction is independent of the data-dimension of the original problem. From a computational point 
of view there are many advantages to this approach. First, the formulas for constructing the mapping of the original 
problem to a matrix problem are all analytic and easy to evaluate, thus computing the relevant reduction is fast. 
Second, the evolution process only involves matrix multiplications, so many data points can be evolved simultaneously 
and, on a multi-core processor, in parallel. Third the time involved in producing the animations showing how the 
points move in data space scales linearly with the number of dimensions to be displayed. Finally, by introducing an 
m that is different from l/cr^ we allow ourselves the freedom of employing low a, which introduces large numbers of 
minima into V, yet also having a low value for m which guarantees efficient tunneling, thus connecting points that 
may be located in nearby, nearly degenerate potential minima. By using this more general Hamiltonian, we reduce 
the sensitivity of the calculation to the specific choice of a. 

One final point worth making before describing the method of calculation is that the use of Gaussian wave-functions 
to represent data points allows us to develop a number of flexible strategies for handling very large data sets. This 
issue will be addressed below. 



The Calculation Method 



Before discussing how this method works in practice we will give a brief outline of the details of the general procedure. 
We begin by assuming that there are n-data points that we wish to cluster. To these data points we associate n-states, 
Itpi). These states are n Gaussian wave- functions such that the i*'^ Gaussian is centered on the coordinates of the i*'^ 
data point. These states form a basis for the vector space within which we calculate the evolution of our model. 

Let us denote by N, the nx n matrix formed from the scalar products 

= m^j), (10) 

and by H, the n x n-matrix 

Hij = mH\i;j), (11) 

and by Xij the matrix of expectation values 



= {tpi\x\t{jj). 



(12) 
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The calculation process can be described in five steps. First, begin by finding the eigenvectors of the symmetric 
matrix N which correspond to states having eigenvalues larger than some pre-assigned value; e.g., 10~^. These 
vectors are linear combinations of the original Gaussians which form an orthonormal set. Second, compute Ti. in 
this orthonormal basis, H*^ . Do the same for Xi j. Fourth, find the eigenvectors and eigenvalues of Ti*'', construct 
lipiit)) = e~**^ lip), that is the solution to the reduced time dependent Schrodinger problem 

i^^\^,(t))^n'nMt)) (13) 

such that \tpi{t — 0)) = Finally, construct the desired trajectories 

{x,{t)) = (V'.le^^'Xe-'*^'^,) (14) 

by evaluating this expression for a range of t and use them to create an animation. Stop the animation when clustering 
of points is obvious. 

It is clear that restricting attention to the truncated Hamiltonian perforce loses some features of the original 
problem, however its advantage is that we can derive analytic expressions for all operators involved (see Appendices 
A and B). As a result, the numerical computations can be done very quickly. Experience has shown that as far as 
clustering is concerned this approximation causes no difhculties. 



Example: Ripley's Crab data 

To test our method we apply it to a five-dimensional dataset with two-hundred entries, used in Ripley's text book 
[lo| . This dataset records five measurements made on male and female crabs that belong to two different species. 
This dataset has been used in the original paper on quantum clustering by Horn and Gottlieb 0]. It is being used here 
to allow readers to compare the two techniques. Applications to other data sets will be discussed below. Our main 
motivation is to provide a simple example which exhibits the details of the DQC method. In particular, we wish to 
show that the simplest computational scheme for implementing the general program captures the essential features 
of the problem and does as well as one can reasonably expect to do. 

The data is stored in a matrix M which has 200 rows and 5 columns. Following an often- used dimensional reduction 
method, we preprocess our data with a singular-value decomposition 

M = USV^, (15) 

where J7 is a unitary 200 x 200 matrix and S is the 200 x 5 matrix of singular values, the latter occurring on the 
diagonal of its upper 5x5 entries. The sub-matrix of U consisting of the first five columns, the so-called five Principal 
Components (PCs), can be thought of as assigning to each sample a unique point in a five-dimensional vector space. 
We may study the problem in the full five-dimensional space, or within any subspace by selecting appropriate principal 
components. In [5] QC was applied to this problem in a 2-dimensional subspace, consisting of PC2 and PCS. In what 
follows we will discuss the application of DQC to the 3-dimensional data composed of the first three PCs (although 
there would be no change in the results if we used all five dimensions). In order to give the reader some feeling 
for how the quantum potential associated with the data looks in a simple case, we have included Fig[l] where we 
exhibit the original data points (different colors of points correspond to known classes), placed upon the associated 
two-dimensional quantum potential. As is obvious from the plot the minima of the potential function do a very 
good job of capturing the different clusters. Moreover, letting data points roll downhill to the nearest minimum will 
produce a reasonable separation of the clusters. 

Clearly, when we restrict attention to the first three PCs, the rows of the matrix obtained by restricting U to its 
first three columns are not guaranteed to be normalized to unity. Hence we employ the conventional approach of 
projecting all points onto the unit sphere [l3]. 

In what follows we study the temporal behavior of the curves {xi{t)), for all i. Henceforth we will refer to this 
as the "motion of points" . Figure [2] shows the distribution of the original data points plotted on the unit sphere in 
three dimensions. This is the configuration before we begin the dynamic quantum evolution. To visually display the 
quality of the separation we have colored the data according to its known four classes, however this information is 
not incorporated into our unsupervised method. To begin with, we see that the two species of crabs ((red, blue) and 
(orange, green)) are fairly well separated; however, separating the sexes in each species is problematic. The middle 
plot in Figure [2] shows the distribution of the points after a single stage of quantum evolution, stopped at a time when 



FIG. 1: This is a plot of the quantum potential function for the two-dimensional problem of the Ripley's Crab Data where the 
coordinates of the data points are chosen to be given by the second and third principal components. The four known classes 
of data points are shown in different colors and are placed upon the potential surface at their original locations. 




FIG. 2: The left hand plot shows three-dimensional distribution of the original data points before quantum evolution. The 
middle plot shows the same distribution after quantum evolution. The right hand plot shows the results of an additional 
iteration of DQC. The values of parameters used to construct the Hamiltonian and evolution operator are: g = 0.07 and 
m — 0.2. Colors indicate the expert classification of data into four classes, unknown to the clustering algorithm. Note, small 
modifications of the parameters lead to the same results. 

points first cross one another and some convergence into clusters has occurred. It is immediately obvious that the 
quantum evolution has enhanced the clustering and made it trivial to separate clusters by eye. Once separation is 
accomplished, extracting the clusters can be performed by eye from the plots or by any conventional technique, e.g. 
k- means. 

An alternative way of displaying convergence is shown in Figure [31 where we plot the Euclidean distance from the 
first point in the dataset to each of the other points. The clusters lie in bands which have approximately the same 
distance from the first point. It is difficult to get very tight clusters since the points, while moving toward cluster 
centers, oscillate around them, and arrive at the minima at slightly different times. Given this intuition, it is clear 
that one way to tighten up the pattern is to stop DQC evolution at a point where the clusters become distinct, and 
then restart it with the new configuration, but with the points redefined at rest. We refer to this as iterating the DQC 
evolution. The right-hand plots in Figure [2] and Figure [3] show what happens when we do this. The second stage of 
evolution clearly tightens up the clusters significantly, as was expected. 

By the end of the second iteration, there can be no question that it is a simple matter to extract the clusters. As 
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FIG. 3: A plot of Euclidean distance of each point i from the first data point. Again, the left hand plot shows the distances 
for the initial distribution of points. The middle plot shows the same distances after quantum evolution. The right-hand plot 
shows results after another iteration of DQC. The numbering of the data-points is ordered according to the expert classification 
of these points into four classes containing 50 instances each. 



is quite evident, clustering does not agree completely with the expert classification, i.e. points with different colors 
may be grouped together. This is, however, the best one can do by color-blind treatment of the information provided 
in the data-matrix. As we already noted the full 5-dimensional study of the crab data-set can proceed in the same 
manner, although it does not lead to new insights. 



Dynamic Distances 



The fact that data-points of different classes happen to lie close to each other in the data-matrix can be due to 
various factors: errors in data measurements, errors in the expert assignment to classes, true proximity of data- 
points in spite of differences of origin (extreme example would be similarities of phenotypes in spite of differences in 
genotypes) or - the simplest possibility - the absence of some discriminative features in the feature-space that spans 
the data measurements. But there is another important conceptual message to be learned here - clustering and/or 
classification may not capture all the interestinglessons that may be derived from the data. A similar message is 
included in the Diffusion Geometry approach 0,13] that advocates measuring diffusion-distances among points rather 
than Euclidean ones. Diffusion distances are influenced by the existence of all other points. In our DQC analysis this 
may be replaced in a straightforward manner by defining dynamic distances among points 

d,,,(0 = ||(f.W>-(f,(t))|| (16) 

with the norm being Euclidean or any other suitable choice. 

Clearly di_j{0) is the geometric distance as given by the original data-matrix or by its reduced form that is being 
investigated. As DQC evolves with time rfi.j(i) changes, and when some semi-perfect clustering is obtained, it will 
be close to zero for points that belong to the same cluster. Figure [3] shows this change in time for all di,i{t) in the 
crab-data example studied above. It is quite obvious that, in addition to the few cases in which clustering disagrees 
with classification, there are many intermediate steps where different data-points are close to each other in spite of 
eventually evolving into different clusters and belonging to different classes. Thus a close scrutiny of the dynamic 
distances matrix dij (t) may lead to interesting observations regarding the relationships among individual pairs of 
points in the original data, a relationship that is brought out by DQC as result of the existing information about all 
other data-points. It may be used to further investigate the reason for such proximities, along any one of the lines 
mentioned above, and thus may lead to novel insights regarding the problem at hand. 



Analysis of Large Data Sets 

There are many scientific and commercial fields, such as cosmology, epidemiology, and risk- management, where 
the data sets of interest contain many points, often also in large numbers of dimensions. We have already discussed 
how to overcome the problem of large dimensions. Dealing with large number of points requires yet a new angle. In 
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problems of this sort it is clear from the outset that, especially on a PC, diagonalizing matrices which are larger than 
2000 X 2000 is computationally intensive. It is obvious that using brute force methods to evolve sets of data having 
tens of thousands of points simply won't work. The solution to the problem of dealing with sets of data containing 
tens of thousands of entries each with N features, lies in the fact that the SVD decomposition maps the data into an 
iV-dimensional cube, and the fact that the data points are represented by states in Hilbert space rather than A^-tuples 
of real numbers. Since there is a literature on ways to do SVD decomposition for large sets of data, we will not address 
this point here. What we do wish to discuss is how to exploit the representation of data points as states in Hilbert 
space in order to evolve large sets of data. 

The trick is to observe that since Gaussian wavcfunctions whose centers lie within a given cube have non- vanishing 
overlaps, as one chooses more and more wavcfunctions one eventually arrives at a situation where the states become 
what we will refer to as essentially linearly dependent. In other words, we arrive at a stage at which any new 
wavefunction added to the set can, to some predetermined accuracy, be expressed as a linear combination of the 
wavcfunctions we already have. Of course, since quantum mechanical time evolution is a linear process, this means 
that this additional state can be evolved by expressing it as a linear combination of the previously selected states and 
evolving them. Since computing the overlap of two Gaussians is done analytically(Appendix B) determining which 
points determine the set of maximally essentially linearly independent states for the problem is easy. Typically, even 
for data sets with 35, 000 points, this is of the order of 1000 points. This works because, as we have already noted, we 
don't need high accuracy for DQC evolution. The quality of the clustering degrades very slowly with loss in accuracy. 
Thus, we can compute the time evolution operator in terms of a well chosen subset of the data and then apply it to 
the whole set of points. [31 

To demonstrate this versatility of DQS we have analyzed a set of 35,213 points in 20 dimensions [13] This data- 
set definitely shows some non-trivial variations in density that can be made apparent by a visual inspection of the 
data plotted in different dimensional combinations. However, DQC is needed to obtain good visual deciphering of 
the different structures. Selecting a sub-set of 1200 points, whose Gaussians we consider to be a set of essentially 
linearly independent states for the problem, we construct Ti*''. By expanding the remaining states in terms of these 
1200 states we can easily evaluate the DQC evolution of all 35,213 points. The results are displayed in Figure 4. It 
seems very clear how the structures develop with DQC. Using the last DQC stage it is possible to identify the main 
structures, and assign each substructure a different color. One can then examine the colored version of a plot of 
the individual data points, discern the structures that belong together, and follow the DQC devlopment tracing out 
dynamics distances between the different points and structures in all dimensions. 



Data exploration involves not only the instances, or data-points, but also the features (coordinates) with which the 
instances are defined. By performing SVD, and selecting a sub-set of coordinates, we define superpositions of the 
original features within which we search for clustering of the instances. In problems with very many features, it is 
adventageous to also perform some feature filtering, employing a judicious selection of subsets of the original features. 
Clearly, the effectiveness of preprocessing data using some method for selecting important features is well appreciated. 
What we wish to show in this discussion is how easily one distinguishes the effects of feature filtering in our visual 
approach and how easy it is, in problems where one has an expert classification, to see if the unsupervised method 
used to select important features is working well. Furthermore, we wish to show the power of combining iterations 
of an SVD based feature filtering algorithm in conjunction with iterations of DQC. To do this we will show what 
happens when one applies these ideas to the dataset of Golub et al. (4|] . 

The Golub et.al. dataset contains gene chip measurements on cells from 72 leukemia patients with two different 
types of Leukemia, ALL and AML. The expert identification of the classes in this data set is based upon dividing the 
ALL set into two subsets corresponding to T-cell and B-cell Leukemia. The AML set is divided into patients who 
underwent treatment and those who did not. In total the Affymetrix GeneChip used in this experiment measured 
the expression of 7129 genes. The feature filtering method we employ is based on SVD-entropy, and is a simple 
modification of a method introduced by Varshavsky et al.^3\ and applied to the same data. 




The method begins by computing the SVD-based entropy 1] of a dataset M ( matrix of n instances by m features of 
Eq. 15) based on the eigenvalues Sj of its diagonal matrix S. Defining normalized relative variance values Vj — y^^-r, 



Interplay of Feature Selection with DQC 
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FIG. 4: A plot of the first three principal components for a large data-set, comprising 35,213 points in 20 dimensions, before 
and after DQC evolution. The potential was determined from the full data-set and evolved using a sub-set of 1200 points, 
whose Gaussians serve as an essentially linearly set of independent states. Three stages of DQC development are shown. The 
coloring was decided upon by selecting the most obvious clusters from the evolved data and assigning colors to them. The dark 
blue points correspond to points that we did not bother to assign to clusters. The purpose of coloring is to be able to look 
at the points in the original data, discern those that belong to common structures, and follow their dynamic distances under 
DQC evolution. 



the dataset entropy is defined through 

where r is the rank of the data-matrix, typically much smaller than m . Given the dataset entropy of the matrix ill, 
define the contribution of the i*'' feature to the entropy using a leave-one-out comparison; i.e., for each feature we 
construct the quantity 

CEi = i?(M(„xm)) - -E^(^(nx(m-l))) (18) 

where the second entropy is computed for the matrix with the i*'' feature removed. Our filtering technique will be to 
remove all features for which CEi < 0. 

Figure O displays the raw data in the 3-dimensional space defined by PCs 2 to 4, and the effect that DQC has on 
these data. In Figure [H we see the result of applying feature filtering to the original data, represented in the same 
3-dimensions, followed by DQC evolution. Applying a single stage of filtering has a dramatic effect upon clustering, 
even before DQC evolution. The latter helps sharpening the cluster separation. 

Figure[7]shows the results of three iterations of SVD-entropy, before and after DQC evolution. These plots, especially 
the after DQC pictures, show dramatic clustering, especially for the blue points. With each stage of filtering we see 
that the blue points cluster better and better, in that the single red outlier separates from the cluster and the cluster 
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FIG. 5: The left hand picture is the raw data from the Affymetrix Chip plotted for principal components 2,3,4. Clearly, without 
the coloring it would be hard to identify clusters. The right hand picture is the same data after DQC evolution using a = 0.2 
and a mass m = 0.01. The different classes are shown as blue, red, green and orange. 




FIG. 6: The left hand plot is the Golub data after one stage of SVD-entropy based filtering, but before DQC evolution. The 
right hand plot is the same data after DQC evolution. 

separates more and more from the other points. The blue points are what we will refer to as an obviously robust cluster 
which has been identified in early stages of filtering. If one continues iterating past the fifth stage, however, the clear 
separation of the blue points from the others begins to diminish. Thus we see that the SVD-entropy based filtering, 
in trying to enhance the clumping of the red points, starts throwing away those features which make the blue cluster 
distinct. Since this effect is quite pronounced we would say that features that are important to distinguishing the 
blue cluster from the others begin to be removed at the sixth and higher iterations of filtering. This is, of course, just 
what we are looking for, a way of identifying those features which are important to the existing biological clustering. 
Out of the original 7129 features, we have reduced ourselves to 2766 features by the fifth iteration. In going from 
step five to step six this gets further reduced to 2488 features, so we could begin searching among the 278 eliminated 
features to isolate those most responsible for the separation of the blue cluster from the others. Instead, we will take 
another track and, since it is so robust and easily identified, remove the blue cluster from the original data and repeat 
the same process without this cluster. The idea here is that now the SVD-entropy based filtering will not be pulled 
by the blue cluster and so it will do a better job of sorting out the red, green and orange clusters. As we will see, this 
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FIG. 7: The left hand plot is the data after three stages of SVD-entropy based filtering, but before DQC evolution. The right 
hand plot is the same data after DQC evolution. 



is in fact the case. 




FIG. 8: The left hand plot is what the starting data looks like if one first removes the blue points and does one stage of 
SVD-entropy based filtering. The right hand plot is what the starting data looks like after three stages of filtering. 

In Figure [8] we see a plot of what the starting configurations look like if one takes the original data, removes the 
identified blue cluster and re-sorts the reduced data set according to the SVD-entropy based filtering rules. The left 
hand plot is what happens if one filters a single time, removing those features, i, whose one-left-out comparison, CEi, 
is less than or equal to zero. The right hand plot shows what happens if one repeats this procedure two more times, 
each time removing features for which CEi < 0. There is no problem seeing that each iteration of the SVD-entropy 
based filtering step improves the separation of the starting clusters. By the time we have done five SVD-entropy 
based filtering steps the red, green and orange clusters are distinct, if not obviously separated. 

Finally, to complete our discussion, we show Figure [9l This figure shows the results of doing five iterations of 
the SVD-entropy based filtering and following that with three stages of DQC evolution. The dramatic clustering 
accomplished by DQC evolution makes it easy to extract clusters. Note however, that in the second plot we see what 
we have seen throughout, that the red points first form two distinct sub-clusters which only merge after two more 
stages of DQC evolution. This constant repetition of the same phenomenon, which is only made more apparent by 
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FIG. 9: The left hand plot is what the starting data looks like if one first removes the blue points and does five stages of 
SVD-entropy based filtering. The right hand plot is what happens after one stage of DQC evolution. The bottom plot is the 
final result after iterating the DQC evolution step two more times. At this point the clusters are trivially extracted. 

SVD-entropy based filtering, is certainly a real feature of the data. It presumably says that what appears to be a 
sample of a single type of cell at the biological level is in reality two somewhat different types of cells when one looks 
at gene expression. A measure of the success of clustering is given by the Jaccard score [l3| which, for this result is 
0.762, higher than the value 0.707 obtained by [13]. 

Summary and Outlook 

We have proposed a dynamical method for exploring proximity relationships among data-points in large spaces. 
Starting with the potential function of quantum clustering Q we have shown how to embed it into a dynamical theory 
so as to provide a visual exploratory tool. Formulating the theoretical treatment using coherent (Gaussian) states 
allows us to derive analytic expressions for all necessary calculations of the temporal evolution. This allows us to 
treat quite complicated data and put them into a visual framework that can be easily manipulated by the user who 
wishes to search for structures in the data. We have tested the system on random data to make sure that it does not 
produce unwarranted clustering structures. 

Throughout this paper we represent the DQC evolution of the Gaussians associated with the original data-points, 
by following the centers of the evolving wave-functions. It should be noted that there is more information to be gained 
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from the full wave-function of a data-point: it is expected to expand, at times, and cover a large fraction of the domain 
of the cluster with which this point is associated. It may also tunnel into neighboring clusters with which the point 
has small dynamic distances. We expect this notion to be particularly useful when the data may be better described 
in terms of 'elongated clusters', i.e. when cluster cores are not points but lines (e.g. a ring) or higher- dimensional 
manifolds. Note that our methodology is not confined to potentials that have only well-separated minima. 

We have discussed the virtues of combining DQC with some preprocessing tools. The first was SVD, which was used 
to limit the range of the data values and to allow us to do some dimensional reduction. While dimensional reduction 
is a must for handling data in very large dimensions, and it helps to remove noise from the data, we wish to point 
out that DQC can handle a large number of features without much difficulty. The computational complexity of the 
problem is controlled by the number of data-points, since this defines the size of the matrix to be exponentiated. The 
computational cost associated with keeping more features is only related to computing the matrices associated with 
multiplying a wave-function by a given coordinate. This is a one time cost. The computational cost of computing 
the values of these operators only grows linearly with the number of features. Clearly it is possible to avoid those 
costs by keeping a large number of features when constructing the quantum potential, V{x), and plotting a much 
smaller number of features when constructing the animations. Experience has shown that after one stage of DQC 
evolution, clustering which occurs because of structures in V{x) that are only seen in features that are not plotted 
in the animations becomes readily visible in those plots that we do construct. This aspect of DQC allows us to 
avoid some of the problems associated with using SVD to strongly reduce the number of dimensions. In addition to 
dimensional reduction based upon simply making an SVD decomposition of the data, we discussed one scheme for 
selecting individual features that are judged to be relevant to the data at hand. Since our problem is unsupervised, we 
employed a feature filtering method that depends on the contribution of the features to SVD-entropy. The examples 
showed that the visual nature of DQC made it easy to judge the effectiveness of feature filtering, especially after 
iterative applications of DQC evolution. 

We have already noted, that for sets of data containing entries with a very large number of features, DQC has 
the computational advantage that once one has formed the Hamiltonian of the system, the computational problem is 
carried out using a matrix which has no more rows and columns than the number of data points. Moreover, we have 
seen that the simplest reduction of the analytic problem of assigning data points to minima of the multi-dimensional 
potential function works remarkably well. Going beyond the truncation procedure explained in Appendix B, while 
easily doable, seems unnecessary for most problems, and this allows us to greatly speed up the computations. In 
our analysis we went on to discuss the case of data sets containing large numbers of points. It turns out that, using 
our Hilbert space representation of data-points, we can naturally select a small set of points whose Gaussians span 
efficiently the entire Hilbert space. These Gaussian are then used as the basis for calculating the DQC evolvement of 
all points. It is quite obvious from the example displayed in Fig. 4 how well these properties of DQC can be employed 
to discern structures in the large data-set under consideration. 

Finally, we wish to observe that the DQC methods described in this paper can be easily extended to general 
classification problems that are usually resolved by supervised machine learning methods. The point is that given a 
training set, i.e., a data set that has been fully resolved by DQC once the appropriate stages of dimensional reduction 
and feature filtering has been applied, then one can use this set to classify new data. There are two different ways one 
can accomplish this task. In the first approach we use the fact that the training set has been successfully clustered 
to assign distinct colors to points that lie in the training set, so that they may be visually identify in all subsequent 
studies. Once this has been done, the classification of new data points can been accomplished in two steps. First, 
reduce the SVD matrix containing both the training set and the new data points (using the previously determined 
features) to an appropriate dimension, and construct the QC potential for the full data set including the training 
set. Next, apply DQC to study the evolution of the full system using the new QC potential and see how the new 
points associate themselves with the points in the training set. Note, as always, both the intermediate dynamics and 
eventual coalescence of the full set into clusters can give useful information about the full data set. The fact that 
the old points have been colored according to the original classification scheme makes it possible to see if the SVD 
reduction of the full data set (training set plus new data) distorts the original classification. If this happens, i.e. if 
the original points fail to cluster properly, then one can go back and use the tools of feature filtering, etc. to analyze 
what has changed. This sort of visual identification of aspects of the data which distort clustering was already used 
in the case of the leukemia data set to see that the existence of a strong cluster can distort the clustering of the 
remaining data. Once this easily identified cluster was removed from the data set the clustering of the remaining data 
was significantly improved. 

The second approach, which is necessary if the dataset contains many entries and the training set is itself large, is 
to use only the training set to generate the quantum potential and the exponential of the Hamiltonian. Next, as we 
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already discussed, use this operator to evolve the full dataset, including the training set. In this case the training set 
is guaranteed to cluster as before and we can categorize the new points according to how they associate with known 
clusters in the training data. 



APPENDIX A. USEFUL OPERATOR IDENTITES 



Using conventional quantum-mechanical notation we represent the Gaussian wave function by 

|^) = (V^a)~^e--'/2a^, (19) 



where we adopted Dirac's bra and ket notation [8| to denote = ipi^) ^-nd ("01 = i'i^Y ■ Employing the operators 

_d_ 

dx 



X and P = ^-j- obeying the commutation relations [x,p\ = i, we define the annihilation operator 



obeying Ao-lcr) = 0. Its Hcrmitian adjoint creation operator A\ = ~* -P + ^ obeys [Aq., Ajj.] = 1. 

We will need a few identities to derive the matrix elements we have to calculate. First we note the normal ordering 
identity (meaning rewriting by using the operator commutation relations so that Aj s appear to the right of all s): 

^a{Al+A^) ^ ^aA^ ^aA„ ^21) 

which may be proven by differentiation with respect to a. Next we note that 

e.MAt ^_,(„)^t ^ Y^gio^^^^ j^t , ^ ^ ^ [AlA^]]] . . .]„ =A„~ 9{a) (22) 

n 

which is easily derived by differentiating with respect to g and noting that only the first commutator is non-zero. A 
similar calculation proves the equally useful result: 

^^(Al-A^) ^ g-aV2 gaAt ^-aA^ ^23) 

Now, because the Parzen window estimator is constructed using Gaussian wavefunctions centered about points 
other than a; = 0, it is convenient to have an operator expression which relates the Gaussian centered about a; = to 
the Gaussian centered about x = x. 

Theorem: |ct, x) = e~'^^^ \a) is a normalized Gaussian wave-function centered at a: = x; i.e. 

1 (x-i)2 

\<j,x) = {^/ira)-^ e 5?^. (24) 

This state is known as a coherent state obeying 

Aa-\a; x) — x\a, x) . (25) 

The generalization to Gaussians in any number of dimensions is straightforward, since they are just products of 
Gaussians defined in each one of the different dimensions. 



APPENDIX B. MATRIX ELEMENTS 



The states we start out with \(7,Xi) have norm one and are, in general, linearly independent; however, they are 
not orthogonal to one another. In what follows we will need an explicit formula for the scalar product of any such 
Gaussian \<T,Xi) with another \cr,Xj). This is easily derived given the operator form for the shifted Gaussian derived 
in Appendix A. Thus we find that 



(26) 
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which is needed for computing the matrix of scalar products Nij = (cr, Xi\a, Xj). Similarly, by employing e*f^ x e ^pv = 
X + y we find that 

y| a; \a, x) = e-(^-^)'/4'^' . (27) 

It is straightforward to generalize this derivation to obtain 

{a,y\V{x)\a,x) = e~^^~y^'/''^\a\V{x + \a), (28) 

for any function V{x). Note that this expectation value can be evaluated by expanding y in a Taylor series about 
the point {x + y)/2. The leading term is simply e~'^^^y^ /^'^ V (^^^) and the remaining terms, involving {a\ a;" \a) 
can be evaluated from the identity 

(-i«"^i-)-E7^H-"i-)-E^F^- (29) 

n=0 ' p=0 ^' 

To speed up computations we chose to approximate all expectation values of V{x) by V{^^^), the first term in this 
series. One could obviously get a more accurate approximation to the original problem by including additional terms 
but explicit computation has shown that, for our purposes, this level of accuracy is sufficient. 

The final formula we need to derive is that for 

{a^p" \a,x) = (ab2e-^(^-^) \a) ^ ^I-^ ^-^^^-yf (30) 

With these preliminaries behind us it only remains to describe the mechanics of the DQC evolution process, where 
we evaluate the Hamiltonian truncated to an n x n matrix in the non-orthonormal basis of shifted Gaussians: 

Hi^j = {a,Xi\H\ij,Xj). (31) 

The time evolution of our original states is computed by applying the exponential of the truncated Hamiltonian to 
the state in question; i.e., |cr, x){t) = e~'^*|cr, x). Computing the exponential of the truncated operator is quite simple, 
except for one subtlety: we have defined 7i by its matrix elements between a non-orthonormal set of states. Hence, to 
perform the exponentiation, we first find the eigenvectors and eigenvalues of the metric Nij and use them to compute 
the matrix N^^^'^ . [l^ Then we construct the transformed Ti by 

=E^M.'^'^M^r/^'- (32) 

k,l 

Now we can construct the exponential of this operator by simply finding its eigenvectors and eigenvalues. In order to 
compute the time evolution of one of the original states we simply write them in terms of the orthonormal basis. 

The only step which remains is to explain how we compute the expectation values of the operator x as functions of 
time: we first construct, for each component, the operator 



Xi^j ^ {a,Xi\x\a,Xj) (33) 



and use N^^ to put this into the same basis in which we exponentiate 7i; i.e., construct 



-1/2 

V)KJ yv.v\j LiiiD iiiLU tJ-it; ooiiiic: v.JC\> 

^^.-E^^'^'^M^^^:/^'- (34) 



k,l 
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